{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# isolated_atom calculation style\n",
    "\n",
    "**Lucas M. Hale**, [lucas.hale@nist.gov](mailto:lucas.hale@nist.gov?Subject=ipr-demo), *Materials Science and Engineering Division, NIST*.\n",
    "\n",
    "## Introduction\n",
    "\n",
    "The isolated_atom calculation style evaluates the base energies of all atomic models associated with an interatomic potential. \n",
    "For some potentials, the isolated energy values are necessary to properly compute the cohesive energy of crystal structures.  This also provides a simple test whether a potential implementation is compatible with a version of LAMMPS. \n",
    "\n",
    "### Version notes\n",
    "\n",
    "- 2020-09-22: Notebook first added.\n",
    "\n",
    "### Additional dependencies\n",
    "\n",
    "### Disclaimers\n",
    "\n",
    "- [NIST disclaimers](http://www.nist.gov/public_affairs/disclaimer.cfm)\n",
    "- Some potentials have two cutoffs with atomic energies outside the first being the \"isolated\" energy while outside the second have zero energy.  The first isolated energy values for those potentials can be found using the diatom_scan calculation instead.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Method and Theory\n",
    "\n",
    "The calculation loops over all symbol models of the potential and creates a system with a single particle inside a system with non-periodic boundary conditions.  The potential energy of each unique isolated atom is evaluated without relaxation/integration.\n",
    "\n",
    "The cohesive energy, $E_{coh}$, of a crystal structure is given as the per-atom potential energy of the crystal structure at equilibrium $E_{crystal}/N$ relative to the potential energy of the same atoms infinitely far apart, $E_i^{\\infty}$\n",
    "\n",
    "$$ E_{coh} = \\frac{E_{crystal} - \\sum{N_i E_{i}^{\\infty}}}{N},$$\n",
    "\n",
    "Where the $N_i$ values are the number of each species $i$ and $\\sum{N_i} = N$.\n",
    "\n",
    "For most potentials, $E_i^{\\infty}=0$ meaning that the measured potential energy directly corresponds to the cohesive energy.  However, this is not the case for all potentials as some have offsets either due to model artifacts or because it allowed for a better fitted model.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Demonstration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Setup"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1.1. Library imports\n",
    "\n",
    "Import libraries needed by the Notebook. The external libraries used are:\n",
    "\n",
    "- [numpy](http://www.numpy.org/)\n",
    "\n",
    "- [atomman](https://github.com/usnistgov/atomman)\n",
    "\n",
    "- [iprPy](https://github.com/usnistgov/iprPy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Notebook last executed on 2020-09-22 using iprPy version 0.10.2\n"
     ]
    }
   ],
   "source": [
    "# Standard library imports\n",
    "import os\n",
    "from pathlib import Path\n",
    "import datetime\n",
    "\n",
    "# http://www.numpy.org/\n",
    "import numpy as np\n",
    "\n",
    "# https://github.com/usnistgov/atomman \n",
    "import atomman as am\n",
    "import atomman.lammps as lmp\n",
    "import atomman.unitconvert as uc\n",
    "\n",
    "# https://github.com/usnistgov/iprPy\n",
    "import iprPy\n",
    "\n",
    "print('Notebook last executed on', datetime.date.today(), 'using iprPy version', iprPy.__version__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1.2. Default calculation setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Specify calculation style\n",
    "calc_style = 'isolated_atom'\n",
    "\n",
    "# If workingdir is already set, then do nothing (already in correct folder)\n",
    "try:\n",
    "    workingdir = workingdir\n",
    "\n",
    "# Change to workingdir if not already there\n",
    "except:\n",
    "    workingdir = Path('calculationfiles', calc_style)\n",
    "    if not workingdir.is_dir():\n",
    "        workingdir.mkdir(parents=True)\n",
    "    os.chdir(workingdir)\n",
    "    \n",
    "# Initialize connection to library\n",
    "library = iprPy.Library(load=['lammps_potentials'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Assign values for the calculation's run parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.1. Specify system-specific paths\n",
    "\n",
    "- __lammps_command__ is the LAMMPS command to use (required).\n",
    "\n",
    "- __mpi_command__ MPI command for running LAMMPS in parallel. A value of None will run simulations serially."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "lammps_command = 'lmp_serial'\n",
    "mpi_command = None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.2. Load interatomic potential\n",
    "\n",
    "- __potential_name__ gives the name of the potential_LAMMPS reference record in the iprPy library to use for the calculation.  \n",
    "\n",
    "- __potential__ is an atomman.lammps.Potential object (required)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "potential_name = '1999--Mishin-Y--Ni--LAMMPS--ipr1'\n",
    "\n",
    "# Retrieve potential and parameter file(s)\n",
    "potential = library.get_lammps_potential(id=potential_name, getfiles=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Define calculation function(s) and generate template LAMMPS script(s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.1 run0.template"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('run0.template', 'w') as f:\n",
    "    f.write(\"\"\"#LAMMPS input script that evaluates a system's energy without relaxing\n",
    "\n",
    "<atomman_system_pair_info>\n",
    "\n",
    "thermo_style custom step pe\n",
    "thermo_modify format float %.13e\n",
    "\n",
    "run 0\"\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.2 isolated_atom()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def isolated_atom(lammps_command, potential, mpi_command=None):\n",
    "    \"\"\"\n",
    "    Evaluates the isolated atom energy for each elemental model of a potential.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    lammps_command :str\n",
    "        Command for running LAMMPS.\n",
    "    potential : atomman.lammps.Potential\n",
    "        The LAMMPS implemented potential to use.\n",
    "    mpi_command : str, optional\n",
    "        The MPI command for running LAMMPS in parallel.  If not given, LAMMPS\n",
    "        will run serially.\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    dict\n",
    "        Dictionary of results consisting of keys:\n",
    "        \n",
    "        - **'energy'** (*dict*) - The computed potential energies for each\n",
    "          symbol.\n",
    "    \"\"\"\n",
    "    # Build filedict if function was called from iprPy\n",
    "    try:\n",
    "        assert __name__ == pkg_name\n",
    "        calc = iprPy.load_calculation(calculation_style)\n",
    "        filedict = calc.filedict\n",
    "    except:\n",
    "        filedict = {}\n",
    " \n",
    "    # Initialize dictionary\n",
    "    energydict = {}\n",
    "    \n",
    "    # Initialize single atom system \n",
    "    box = am.Box.cubic(a=1)\n",
    "    atoms = am.Atoms(atype=1, pos=[[0.5, 0.5, 0.5]])\n",
    "    system = am.System(atoms=atoms, box=box, pbc=[False, False, False])\n",
    "\n",
    "    # Get lammps units\n",
    "    lammps_units = lmp.style.unit(potential.units)\n",
    "\n",
    "    # Define lammps variables\n",
    "    lammps_variables = {}\n",
    "\n",
    "    # Loop over symbols\n",
    "    for symbol in potential.symbols:\n",
    "        system.symbols = symbol\n",
    "\n",
    "        # Add charges if required\n",
    "        if potential.atom_style == 'charge':\n",
    "            system.atoms.prop_atype('charge', potential.charges(system.symbols))\n",
    "\n",
    "        # Save configuration\n",
    "        system_info = system.dump('atom_data', f='isolated.dat',\n",
    "                                  potential=potential,\n",
    "                                  return_pair_info=True)\n",
    "        lammps_variables['atomman_system_pair_info'] = system_info\n",
    "        \n",
    "        # Write lammps input script\n",
    "        template_file = 'run0.template'\n",
    "        lammps_script = 'run0.in'\n",
    "        template = iprPy.tools.read_calc_file(template_file, filedict)\n",
    "        with open(lammps_script, 'w') as f:\n",
    "            f.write(iprPy.tools.filltemplate(template, lammps_variables,\n",
    "                                             '<', '>'))\n",
    "        \n",
    "        # Run lammps and extract data\n",
    "        output = lmp.run(lammps_command, lammps_script, mpi_command)\n",
    "        energy = output.simulations[0]['thermo'].PotEng.values[-1]\n",
    "        energydict[symbol] = uc.set_in_units(energy, lammps_units['energy'])\n",
    "    \n",
    "    # Collect results\n",
    "    results_dict = {}\n",
    "    results_dict['energy'] = energydict\n",
    "    \n",
    "    return results_dict"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4. Run calculation function(s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "results_dict = isolated_atom(lammps_command, potential, mpi_command=mpi_command)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dict_keys(['energy'])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results_dict.keys()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5. Report results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5.1. Define units for outputting values\n",
    "\n",
    "- __energy_unit__ is the unit of energy to display values in."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "energy_unit = 'eV'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5.2 Display isolated energies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Ni -3.0970029925997e-11 eV\n"
     ]
    }
   ],
   "source": [
    "for symbol, energy in results_dict['energy'].items():\n",
    "    print(symbol, uc.get_in_units(energy, energy_unit), energy_unit)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
